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Abstract — The sparse signal recovery in tiie standard com- 
pressed sensing (CS) problem requires that the sensing matrix 
be linown a priori. Such an ideal assumption may not be met in 
practical applications where various errors and fluctuations exist 
in the sensing instruments. This paper considers the problem 
of compressed sensing subject to a structured perturbation in 
the sensing matrix. Under mild conditions, it is shown that 
a sparse signal can be recovered by H-i minimization and the 
recovery error is at most proportional to the measurement 
noise level, which is similar to the standard CS result. In the 
special noise free case, the recovery is exact provided that the 
signal is sufficiently sparse with respect to the perturbation 
level. The formulated structured sensing matrix perturbation is 
applicable to the direction of arrival estimation problem, so has 
practical relevance. Algorithms are proposed to implement the 
£i minimization problem and numerical simulations are carried 
out to verify the result obtained. 

Index Terms — Compressed sensing, structured matrix pertur- 
bation, stable signal recovery, alternating algorithm, direction of 
arrival estimation. 



I. Introduction 

Compressed sensing (CS) has been a very active research 
area since the pioneering works of Candes et at. |[T], Q and 
Donoho |3]. In CS, a signal x° E M" of length n is called 
fc-sparse if it has at most k nonzero entries, and it is called 
compressible if its entries obey a power law 



(1) 



where 



is the jth largest entry (in absolute value) of x° 

^ l^°l(2) ^ ■ ■ ■ ^ l^°l(n))' 9 > 1 Cq is a constant 



that depends only on q. Let x'' be a vector that keeps the k 
largest entries (in absolute value) of x° with the rest being 
zeros. If x° is compressible, then it can be well approximated 
by the sparse signal x'^ in the sense that 



\x — X 



< C'k-'i+^'^ 



(2) 



12 — 1 

where C'^ is a constant. To obtain the knowledge of x°, CS 
acquires linear measurements of x° as 

^x° 



y 



e, 



(3) 



where $ e jjmxn j-j^g sensing matrix (or linear operator) 
with typically k < m <^ n, y E M™ is the vector of mea- 
surements, and e e M™ denotes the vector of measurement 
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noises with bounded energy, i.e., jjejlj < e for e > 0. Given 
# and e, the task of CS is to recover x° from a significantly 
reduced number of measurements y. Candes et al. ||T], Q 
show that if x° is sparse, then it can be stably recovered 
under mild conditions on $ with the recovery error being 
at most proportional to the measurement noise level e by 
solving an £i minimization problem. Similarly, the largest 
entries (in absolute value) of a compressible signal can be 
stably recovered. More details are presented in Subsection 
II-B In addition to the li minimization, other approaches that 
provide similar guarantees are also reported thereafter, such 
as IHT |5 j and greedy pursuit methods including OMP |j6j, 
StOMP [7 1 and CoSaMP |8l. 

The sensing matrix ^ is assumed known a priori in standard 
CS, which is, however, not always the case in practical 
situations. For example, a matrix perturbation can be caused 
by quantization during implementation. In source separation 
|j9), pO| the sensing matrix (or mixing system) is usually 
unknown and needs to be estimated, and thus estimation 
errors exist. In source localization such as direction of arrival 
(DOA) estimation (|TT), ||T2) and radar imaging |[T3), |[T4|, the 
sensing matrix (overcomplete dictionary) is constructed via 
discretizing one or more continuous parameters, and errors 
exist typically in the sensing matrix since the true source 
locations may not be exactly on a discretized sampling grid. 

There have been recent active studies on the CS problem 
where the sensing matrix is unknown or subject to an unknown 
perturbation. Gleichman and Eldar f\5\ introduce a concept 
named as blind CS where the sensing matrix is assumed 
unknown]^ In order for the measurements y to determine a 
unique sparse solution, three additional constraints on $ are 
studied individually and sufficient conditions are provided to 
guarantee the uniqueness. Herman and Strohmer p6) analyze 
the effect of a general matrix perturbation and show that the 
signal recovery is robust to the perturbation in the sense that 
the recovery error grows linearly with the perturbation level. 
Similar robust recovery results are also reported in fTT) , p8) . 
It is demonstrated in |jT8), |I9| that the signal recovery may 
suffer from a large error under a large perturbation. In addition, 
the existence of recovery error caused by the perturbed sensing 
matrix is independent of the sparsity of the original signal. 
Algorithms have also been proposed to deal with sensing 
matrix perturbations. Zhu et al. IpOl propose a sparse total 



'The CS problem formulation in jis] is a little different from that in jsj. 
In 1 15 1, the signal of interest is assumed to be sparse in a sparsity basis while 
the sparsity basis is absorbed in the sensing matrix $ in our formulation. The 
sparsity basis is assumed unknown in \\5\ that leads to an unknown sensing 
matrix in our formulation. 
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least-squares approach to alleviating the effect of perturbation 
where they explore the structure of the perturbation to improve 
recovery performance. Yang et al. |12| formulate the off-grid 
DOA estimation problem from a sparse Bayesian inference 
perspective and iteratively recover the source signal and the 
matrix perturbation. It is noted that existing algorithmic results 
provide no guarantees on signal recovery accuracy when there 
exist perturbations in the sensing matrix. 

This paper is on the perturbed CS problem. A structured 
matrix perturbation is studied with each column of the per- 
turbation matrix being a (unknown) constant times a (known) 
vector which defines the direction of perturbation. For cer- 
tain structured matrix perturbation, we provide conditions for 
guaranteed signal recovery performance. Our analysis shows 



CS. Sectionjnijpresents the main results of the paper as well as 
some discussions and a practical application in DOA estima- 



that robust stability (see definition in Subsection II-A i can be 
achieved for a sparse signal under similar mild conditions as 
those for standard CS problem by solving an £i minimization 
problem incorporated with the perturbation structure. In the 
special noise free case, the recovery is exact for a sufficiently 
sparse signal with respect to the perturbation level. A similar 
result holds for a compressible signal under an additional 
assumption of small perturbation (depending on the number 
of largest entries to be recovered). A practical application 
problem, the off-grid DOA estimation, is further considered. It 
can be formulated into our proposed signal recovery problem 
subject to the structured sensing matrix perturbation, showing 
the practical relevance of our proposed problem and solution. 
To verify the obtained results, two algorithms for positive- 
valued and general signals respectively are proposed to solve 
the resulting nonconvex £i minimization problem. Numerical 
simulations confirm our robustly stable signal recovery results. 

A common approach in CS to signal recovery is solving 
an optimization problem, e.g., £i minimization. In this con- 
nection, another contribution of this paper is to characterize a 
set of solutions to the optimization problem that can be good 
estimates of the signal to be recovered, which indicates that it 
is not necessary to obtain the optimal solution to the optimiza- 
tion problem. This is helpful to assess the "effectiveness" of an 
algorithm (see definition in Subsection |ni-E| i, for example, the 
£p ip < 1) minimization | 2Tj , |j22j in standard CS, in solving 
the optimization problem since in nonconvex optimization the 
output of an algorithm cannot be guaranteed to be the optimal 
solution. 

Notations used in this paper are as follows. Bold-case letters 
are reserved for vectors and matrices. ||a;||Q denotes the pseudo 
£q norm that counts the number of nonzero entries of a vector 
X. ||a;||j^ and ||a;||2 denote the £i and £2 norms of a vector x 
respectively. ||A||2 and ||A||p are the spectral and Frobenius 
norms of a matrix A respectively, x^ is the transpose of a 
vector X and is for a matrix A. Xj is the jth entry of 
a vector x. is the complementary set of a set T. Unless 
otherwise stated, xt has entries of a vector x on an index 
set T and zero entries on T'^. diag [x) is a diagonal matrix 
with its diagonal entries being entries of a vector a;. is the 
Hadamard (elementwise) product. 

The rest of the paper is organized as follows. Section [III first 
defines formally some terminologies used in this paper and 
then introduces existing results on standard CS and perturbed 



tion. Section IV introduces algorithms for the £1 minimization 
problem in our considered perturbed CS and their analysis. 
Section |V] presents extensive numerical simulations to verify 
our main results and also empirical results of DOA estimation 
to support the theoretical findings. Conclusions are drawn in 



Section VI Finally, some mathematical proofs are provided in 



Appendices. 

II. Preliminary Results 

A. Definitions 

For the purpose of clarification of expression, we define 
formally some terminologies for signal recovery used in this 
paper, including stability in standard CS, robustness and robust 
stability in perturbed CS. 

Definition 1 ( [1]): In standard CS where $ is known a 
priori, consider a recovered signal x of x° from measurements 
y — ^x° + e with ||e||2 < e. We say that x achieves stable 
signal recovery if 

\\x - x°\\^ < Cf Q^tb^ 

holds for compressible signal x° obeying ([TJ and an integer 
k, or if 

Wx-x" 



holds for /c-sparse signal x°, with nonnegative constants , 

Definition 2: In perturbed CS where = A + E with A 
known a priori and E unknown with ||-E|lp < 77, consider a 
recovered signal x of x° from measurements y = ^x" + e 
with ||e||2 < e. We say that x achieves robust signal recovery 
if 

\\X - X°\\^ < C[''*fc-«+l/2 ^ (jrbt^ ^ (jrbt^ 

holds for compressible signal x° obeying ([T]i and an integer 
fc, or if 



\x — X 



rbt^ 
3 



holds for fc-sparse signal x°, with nonnegative constants C[ , 
C^''* and Cl^K 

Definition 3: In perturbed CS where ^ — A + E with A 
known a priori and E unknown with ||-E||p < rj, consider a 
recovered signal x of x° from measurements y = ^x" + e 
with ||e||2 < e. We say that x achieves robustly stable signal 
recovery if 

\\x-x°\\^<Cr (r;)fc-^+i/2 + C2- (r7)e 

holds for compressible signal x° obeying ([T]i and an integer 
fc, or if 



\x — X 



<C;^rj)e 



holds for fc-sparse signal x°, with nonnegative constants C['*, 
C2'' depending on rj. 
Remark 1: 

(1) In the case where x° is compressible, the defined stable, 
robust, or robustly stable signal recovery is in fact for 
its fc largest entries (in absolute value). The first term 
O (fc^9+^/'^) in the error bounds above represents, by (j2ji. 
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the best approximation error (up to a scale) that can be 
achieved when we know everything about x° and select 
its k largest entries. 

(2) The Frobenius norm of E, ||-E||p, can be replaced by 
any other norm in Definitions |2] and |3] since the norms 
are equivalent. 

(3) By robust stability, we mean that the signal recovery is 
stable for any fixed matrix perturbation level t] according 
to Definition [3] 

It should be noted that the stable recovery in standard CS 
and the robustly stable recovery in perturbed CS are exact in 
the noise free, sparse signal case while there is no such a 
guarantee for the robust recovery in perturbed CS. 



B. Stable Signal Recovery of Standard CS 

The task of standard CS is to recover the original signal x° 
via an efficient approach given the sensing matrix acquired 
sample y and upper bound e for the measurement noise. This 
paper focuses on the £i norm minimization approach. The 
restricted isometry property (RIP) (23] has become a dominant 
tool to such analysis, which is defined as follows. 

Definition 4: Define the /c-restricted isometry constant 
(RIC) of a matrix denoted by ($), as the smallest 
number such that 



(1-4 m\\v\\i<\\^v 



l2<(l+4 (*))||^||2 



holds for all /c-sparse vectors v. ^ is said to satisfy the fc-RIP 
with constant Sk (*) if 4 (*) < 1. 

Based on the RIP, the following theorem holds. 

Theorem 1 ( [4]): Assume that (*) < \/2 — 1 and 
||e|l2 — Then an optimal solution x* to the basis pursuit 
denoising (BPDN) problem 



a;||-^ , subject to ||y — ^xW^ < e 



satisfies 



X — X 



I2 < Q*'^fc"l/2 



\x — X 



CI 



(4) 



(5) 



where C^*'' = 



_ 2[l+(x/2-l)^2fc(^)] 



CI 



td 



■"0 ^ l-(y2+l)52fc(*) ' ~ l-(V2+l)52fc(*)' 

Theorem[I]states that a fc-sparse signal x° {x^ = x°) can be 
stably recovered by solving a computationally efficient convex 
optimization problem provided 52k (^) < \/2 — 1. The same 
conclusion holds in the case of compressible signal x" since 



1^ <c;'fc-«+i/2 



(6) 



ng to (TTJ and ^ with C^' being a constant. In the 
noise free, fc-sparse signal case, such a recovery is 



according 

special noise tree, fe-sparse signal case, such a recovery 
exact. The RIP condition in Theorem [T] can be satisfied 
provided m > O (fclog [n/k]) with a large probability if the 
sensing matrix ^ is i.i.d. subgaussian distributed |24|. Note 
that the RIP condition for the stable signal recovery in standard 
CS has been relaxed in p5) , p6) but it is beyond the scope 
of this paper. 



C. Robust Signal Recovery in Perturbed CS 

In standard CS, the sensing matrix $ is assumed to be 
exactly known. Such an ideal assumption is not always the 
case in practice. Consider that the true sensing matrix is 
— A + E where A € jijmxn ^j^g lojown nominal 
sensing matrix and E G M™^" represents the unknown 
matrix perturbation. Unlike the additive noise term e in the 
observation model in (|3]l, a multiplicative "noise" Ex° is 
introduced in perturbed CS and is more difficult to analyze 
since it is correlated with the signal of interest. Denote ||-E||2'^^ 
the largest spectral norm taken over all A; -column submatrices 
of E, and similarly define ||$||2'^\ The following theorem is 
stated in fT6l. 



Theorem 2 ( /[76|/j.- Assume that there exist constants 



e and eE.x° such that 



< £ 



(k) 



e||2 < e and 



\Ex°\ 



Assume that S2k (*) < 



V2 



1+e 



(2k)\ 



- 1 



and ||a;°||Q < k. Then an optimal solution x* to the BPDN 
problem with the nominal sensing matrix A, denoted by N- 
BPDN, 



Ax\\^ <e 



min||a;|L, subject to \\y 

X 

achieves robust signal recovery with 

\\x* -x°\\^ < CP'^e + C^'^E 



(7) 



(8) 



Vi+A-2fc(#)(i+4"':^) 



where CP*'' = - 

l-{V2+l)^{l+S2km){l+e'^''l) -l\ 

Remark 2: 

(1) The relaxation of the inequality constraint in ^ from 
e to e + ce.x" is to ensure that the original signal 
03° is a feasible solution to N-BPDN. Theorem |2] is a 
little different from that in fW\, where the multiplicative 
"noise" Ex° is bounded using e^"*^, 4 (^) and ||$a;°||2 
rather than a constant 3,0. 

(2) Theorem |2] is applicable only to the small perturbation 
case where < — 1 since 62k (^) > 0. 

(3) Theorem |2] generalizes Theorem [T] for the fc-sparse signal 
case. As the perturbation E ^ 0, Theorem |2] coincides 
with Theorem [T] for the fc-sparse signal case. 

Theorem |2] states that, for a small matrix perturbation 
E, the signal recovery of N-BPDN that is based on the 
nominal sensing matrix A is robust to the perturbation with the 
recovery error growing at most linearly with the perturbation 
level. Note that, in general, the signal recovery in Theorem 
|2] is unstable according to the definition of stability in this 
paper since the recovery error cannot be bounded within a 
constant (independent of the noise) times the noise level as 



some perturbation occurs. A result on general signals in 1 16| 



is omitted that shows the robust recovery of a compressible 
signal. The same problem is studied and similar results are 
reported in |17J based on the greedy algorithm CoSaMP fSl. 

III. SP-CS: CS Subject to Structured Perturbation 

A. Problem Description 

In this paper we consider a structured perturbation in the 
form E = BA° where B e E™^" is known a priori, A° = 



4 



diag (/3°) is a bounded uncertain term with (3° E [— r, r]" and 
r > 0, i.e., each column of the perturbation is on a known 
direction. In addition, we assume that each column of B has 
unit norm to avoid the scaling problem between B and A° (in 



where 



fact, the D-RIP condition on matrix [A, B] in Subsection III-B 



implies that columns of both A and B have approximately unit 
norms). As a result, the observation model in (|3]l becomes 



y 



BA° 



(9) 



ejlj < e. Given 



with A° = diag(/3°), (3° e [-r,r]" and 
y. A, B, r and e, the task of SP-CS is to recover x° and 
possibly /3° as well. 
Remark 3: 

(1) Without loss of generality, we assume that x, y. A, B 
and e are all in the real domain unless otherwise stated. 

(2) If a;° = for some j e {I,-- - ,n}, then /3° has 
no contributions to the observation y and hence it is 
impossible to recover (3°. As a result, the recovery of /3° 
in this paper refers only to the recovery on the support 
of x°. 



B. Main Results of This Paper 

In this paper, a vector v is called 2fc-duplicately (D-) sparse 
if t; = [t)f , with Vi and V2 being of the same dimension 
and jointly fc-sparse (each being fc-sparse and sharing the same 
support). The concept of duplicate (D-) RIP is defined as 
follows. 

Definition 5: Define the 2fc-duplicate (D-) RIC of a matrix 
denoted by as the smallest number such that 



(1 - hk (*)) ll^^ll' < Il*^ll2 < (1 + ~^2k (*)) 11^^ 



holds for all 2fc-D-sparse vectors v. €> is said to satisfy the 
2fc-D-RIP with constant (*) if hk (*) < 1- 

With respect to the perturbed observation model in (|9]l, let 
^ — [A, B]. The main results of this paper are stated in the 
following theorems. The proof of TheoremJS] is provided in 
Appendix A and proofs of Theorems [4] and [Sjare in Appendix 
B. 

Theorem 3: In the noise free case where e = 0, assume 
that ||3;°||o < k and 54k (*) < 1- Then an optimal solution 
{x* ,f3*) to the perturbed combinatorial optimization problem 



mm 

a;GK",/3G[-r,r] 



IIq , subject to y = (A + SA) a; (10) 
with A = diag {(3) recovers x° and f3°. 



Theorem 4: Assume 



that 



5ik (*) 



< 



x°\\q ^ k and ||e||2 < e. Then 



r^) + l 

an optimal solution {x*,f3*) to the perturbed (P-) BPDN 
problem 

min INIlii subject to \\y — {A + BA.) x\\^ < e 

(11) 

achieves robustly stable signal recovery with 



\x* -a;°||2 < Ce, 
(/3* -/3°)0a;°||2 <Ce 



(12) 
(13) 



C = 



1- (v/2(l + r2) + lU4,.(*) 



[2 + \/rT72 11^112 c] 

Vl - S4k (*) 



Theorem 5: Assume that d^k (*) < (1 + r^) + 1 
and ||e||2 < e. Then an optimal solution {x*,(3*) to the P- 
BPDN problem in ( [TT] ) satisfies that 

\\x* -x^W^ < (Cofc-'/' + Ci) ||a;°-a;'=||^+C2e(14) 
\\if3* ~(3'')Qx% 

< (Cofc-i/'+Ci) ||a;°-a;'=||^+C2e (15) 



where 



Co = 2 



l+(^V2(l + ^2)-l)'^4fe(*) /a, 
Ci = 2V2r54k (*) /a, 
Co = v/l + 11*112^0/5, 



Ci = 



VT+72Ci + 2r 



m2/b 



with a = 1 - (v/2(l + r2) + l) ^4,, (*), b = ^1 - 54k (*) 



and C2 = C, C2 = C with C, C as defined in Theorem [4] 

Remark 4: In general, the robustly stable signal recovery 
cannot be concluded for compressible signals since the error 
bound in ([T4]i may be very large in the case of large per- 
turbation by Ci = 0{r). If the perturbation is small with 
r = O (fc^^/2j^ then the robust stability can be achieved for 
compressible signals by (|6]l provided that the D-RIP condition 
in Theorem |5] is satisfied. 



C. Interpretation of the Main Results 

Theorem [5] states that for a fc-sparse signal x°, it can be 
recovered by solving a combinatorial optimization problem 
provided 54k {^) < 1 when the measurements are exact. 
Meanwhile, /3° can be recovered. Since the combinatorial 
optimization problem is NP-hard and that its solution is 
sensitive to measurement noise |27j, a more reliable approach, 
£1 minimization, is explored in Theorems |4] and |5] 

Theorem |4] states the robustly stable recovery of a fc- 
sparse signal x° in SP-CS with the recovery error being at 
most proportional to the noise level. Such robust stability is 
obtained by solving an £4 minimization problem incorporated 
with the perturbation structure provided that the D-RIC is 
sufficiently small with respect to the perturbation level in terms 
of r. Meanwhile, the perturbation parameter f3° can be stably 
recovered on the support of x". As the D-RIP condition is 
satisfied in Theorem |4] the signal recovery error of perturbed 
CS is constrained by the noise level e, and the influence of the 
perturbation is limited to the coefficient before e. For example, 
if <54fc(*) 0.2, then ||a;* - a;°||2 < 8.48e, 8.50e, ll.Oe 
corresponding to r = 0.01,0.1, 1, respectively. In the special 
noise free case, the recovery is exact. This is similar to that in 
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Standard CS but in contr ast to the existing robust signal recov- an example where x° is A:-sparse. Let z° = 



ery result in Subsection II-C where the recovery error exists 
once a matrix perturbation appears. Another interpretation of 
the D-RIP condition in Theorem |4] i s that the robustly stable 



and 



signal recovery requires that r < ^ ^ yS^k (^P) ^ ij —1 
for a fixed matrix '3/. Using the aforementioned example where 
Sik ('S') — 0.2, the perturbation is required to satisfy r < y/7. 
As a result, our robustly stable signal recovery result of SP- 
CS applies to the case of large perturbation if the D-RIC of 
^ is sufficiently small while the existing result does not as 
demonstrated in Remark |2] 

Theorem [5] considers general signals and is a generalized 
form of Theorem]?] In comparison with Theorem]T]in standard 
CS, one more term Ci ||a;° — a;*^ || ^ appears in the upper 
bound of the recovery error. The robust stability does not hold 
generally for compressible signals as illustrated in Remark]?] 
while it is true under an additional assumption r ^ O (^k^^^^y 
The results in this paper generalize that in standard CS. 
Without accounting for the symbolic difference between 
(^) and 64k ('4'), the conditions in Theorems ]T] an d ]5] 
coincide, as well as the upper bounds in Q and (l4\ for 
the recovery errors, as the perturbation vanishes or equiva- 
lently r — > 0. As mentioned before, the RIP condition for 
guaranteed stable recovery in standard CS has been relaxed. 
Similar techniques may be adopted to possibly relax the D- 
RIP condition in SP-CS. While this paper is focused on the 
£1 minimization approach, it is also possible to modify other 
algorithms in standard CS and apply them to SP-CS to provide 
similar recovery guarantees. 

D. When is the D-RIP satisfied? 

Existing works studying the RIP mainly focus on random 
matrices. In standard CS, ^ has the fc-RIP with constant 5 with 
a large probability provided that m > Csk\og{n/k) and ^ 
has properly scaled i.i.d. subgaussian distributed entries with 
constant Cs depending on 5 and the distribution p4) . The 
D-RIP can be considered as a model-based RIP introduced 



in (28|. Suppose that A, B are mutually independent and 
both are i.i.d. subgaussian distributed (the true sensing matrix 
# = A + BA° is also i.i.d. subgaussian distributed if /3° is 
independent of A and B). The model-based RIP is determined 
by the number of subspaces of the structured sparse signals 
that are referred to as the D-sparse ones in the present paper 
For ^ — [A,B], the number of 2A; -dimensional subspaces for 



2fc-D-sparse signals is ^^J . Consequently, ^ has the 2fc-D- 
RIP with constant 5 with a large probability also provided that 



m > Csklog (n/k) by (j28) Theorem 1] or |29 Theorem 3.3]. 
So, in the case of a high dimensional system and r — > 0, the 
D-RIP condition on \J' in Theorem|4]or|5]can be satisfied when 
the RIP condition on 4> (after proper scaling of its columns) 
in standard CS is met. It means that the perturbation in SP-CS 
gradually strengthens the D-RIP condition for robustly stable 
signal recovery but there exists no gap between SP-CS and 
standard CS in the case of high dimensional systems. 

It is noted that there is another way to stably recover the 
original signal x° in SP-CS. Given the sparse signal case as 



X" 

x° 

it is 2A;-sparse. The observation model can be written as y = 
+ e. Then z° and hence, x°, can be stably recovered 
from the problenj^ 



subject to \\y ■ 



(16) 



provided that (*) < %/2 — 1 by Theorem [l] It looks like 
that we transformed the perturbation into a signal of interest. 



Denote TPS-BPDN the problem in ( 16 1. In a high dimensional 



system, the condition \/2 — 1 requires about twice 

as many as the measurements that makes the D-RIP condition 
5ik ('4') < \/2 — 1 hold by |28 Theorem 1] corresponding 
to the D-RIP condition in Theorem ]4] or ]5] as r 0. As 
a result, for a considerable range of perturbation level, the 
D-RIP condition in Theorem |4] or |5] for P-BPDN is weaker 
than that for TPS-BPDN since it varies slowly for a moderate 
perturbation (as an example, S^^k (*) < 0.414, 0.413, 0.409 
corresponds to r = 0,0.1,0.2 respectively). Numerical simu- 
lations in Subsection |V] can verify our conclusion. 

E. Relaxation of the Optimal Solution 

In Theorem ]5] (Theorem ]4] is a special case), {x*,(3*) is 
required to be an optimal solution to P-BPDN. Naturally, we 
would like to know if the requirement of the optimality is 
necessary for a "good" recovery in the sense that a good 
recovery validates the error bounds in (JT4li and (JT5]l under 



the conditions in Theorem ]5] Generally speaking, the answer 
is negative since, regarding the optimality of (a;*,/3*), only 
< and the feasibility of {x* ,13*) are used in the 

proof of Theorem ]5] in Appendix B. Denote V the feasible 
domain of P-BPDN, i.e.. 



I?-{(a;,/3) :/3e [-r,r]", 

\\V - (A + BA) x\\^ < e with A = diag (/3)}. 



(17) 



We have the following corollary. 

Corollary 1: Under the assumptions in Theorem ]5] any 
{x,/3) E T) that meets ||a;||j^ < ||a;°||]^ satisfies that 



\\x - x°\\^ < (Cok-^l^ + Ci) ||a;° - a; 



2e, 



|(/3-/3°)0a;* 



< 



-1/2 



\x — X 



+ C2e 



with Cj,Cj, j = 0, 1, 2, as defined in Theorem ]5] 

Corollary ]T] generalizes Theorem ]5] and its proof follows 
directly from that of Theorem ]5] It shows that a good recovery 
in SP-CS is not necessarily an optimal solution to P-BPDN. A 
similar result holds in standard CS that generalizes Theorem 
JT] and the proof of Theorem JT] in |4| applies directly to such 
case. 

Corollary 2: Under the assumptions in Theorem ]!] any x 
that meets \\y — AxW^ < e and ||a;||j^ < ||a;°||]^ satisfies that 



\\x - x°\\^ < ||a;° -x'^l 

with C^*^,Cf*'' as defined in Theorem [l] 



(18) 
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A,{e) = a{ 



Fig. 1. Illustration of Corollary|2] The shaded band area refers to the feasible 
domain of BPDN. The triangular area, the intersection of the feasible domain 
and the £i ball {x : \\x\\-^ < ||a;°||-|^ }, is the set of all good recoveries. 



An illustration of Corollary |2] is presented in Fig. [T] 
where the shaded band area refers to the feasible domain 
of BPDN in (j4| and all points in the triangular area, the 
intersection of the feasible domain of BPDN and the £i ball 



{x 



< 



|a;°||j^}, are good candidates for recovery of x°. 



The reason why one seeks for the optimal solution x* is to 
guarantee that the inequality ||a;||j^ < ||a;°||j^ holds since ||a;°||]^ 
is generally unavailable a priori. Corollary |2] can explain why 
a satisfactory recovery can be obtained in practice using some 
algorithm that may not produce an optimal solution to BPDN, 
e.g., rONE-Ll [30l. Corollaries [T]and|2]are useful for checking 
the effectiveness of an algorithm in the case when the output 
cannot be guaranteed to be optimal]^ Namely, an algorithm is 
called effective in solving some £i minimization problem if it 
can produce a feasible solution with its £i norm no larger than 
that of the original signal. Similar ideas have been adopted in 



F. Application to DOA Estimation 

DOA estimation is a classical problem in signal processing 
with many practical applications. Its research has recently been 
advanced owing to the development of CS based methods, e.g., 
£i-SYD ([34). This subsection shows that the proposed SP-CS 
framework is applicable to the DOA estimation problem and 
hence has practical relevance. Consider k narrowband far-field 
sources sj, j — 1, ■ ■ ■ , k, impinging on an m-element uniform 
linear array (ULA) from directions dj with dj £ [0, vr), j = 
I,-- - ,k. Denote 6 = [cos (di) , . . . , cos (dfe)]^ £ (-1,1]''. 
For convenience, we consider the estimation of 9 rather 
than that of d hereafter Moreover, we consider the noise 
free case for simplicity of exposition. Then the observation 
model is y = A (0) s according to p5) , where y G C™ 
denotes the vector of sensor measurements and A {0) e 
denotes the sensing/measurement matrix with respect to 6 with 



into the problem 



^It is hard to incorporate the knowledge /3° G [— 
in {16}. 

It is common when the problem to be solved is nonconvex, such as 
P-BPDN as discussed in Section |IV| and £p (0 < p < 1) minimization 
approaches i21J, [,22] , [31] in standard CS. In addition. Corollaries ^ and [2] 
can be readily extendecito the £p (0 < p < 1) minimization approaches. 



and ai{e,) = ^^cxp {tn {l - ^) 9,}, 
j = 1, • • • , fc, I — 1, • • • , TO, where i = The 
objective of DOA estimation is to estimate 6 given y and 
possibly k as well. Since k is typically small, CS based 
methods have been motivated in recent years. Let 6 = 
{i — 1,| — 1,---,1— be a uniform sampling grid in 
6 range (—1,1] where n denotes the grid number (without 
loss of generality, we assume that n is an even number). In 
existing standard CS based methods 9 actually serves as the set 
of all candidate DOA estimates. As a result, their estimation 
accuracy is limited by the grid density since for some 9j, 
j e {l,-'' jk}, the best estimate of 6j is its nearest grid 
point in 9. It can be easily shown that a lower bound for the 
mean squared estimation error of each 9j is LB = by 
assuming that 9j is uniformly distributed in one or more grid 
intervals. 

An off-grid model has been studied in |T2), that takes 
into account effects of the off-grid DOAs and introduces a 
structured matrix perturbation in the measurement matrix. For 
completeness, we re-derive it using Taylor expansion. Suppose 
9j ^ 9 for some j e {1, • • • 
is the nearest grid point to 



, k} and that 9i^, Ij e {1, • • • , n}, 
'j. By Taylor expansion we have 



ai9,) = a[9i^+b 



(19) 



with b ( 9, 



term with respect to 9j. Denote k = | 

B = K 



a \9ij ,a[9n 
for ( = 1 , • • • ,n 



and Rj being a remainder 
b {9i ] , 



and 



3 ' 

b(9,. 



/3f = K (^9j - 9i^j , xf = sj, ifl^ Ij for any j € {1, • • • , fc} : 
13° = 0, x° = 0, otherwise, 

with Ij e {1, • • • , n} and 9i. being the nearest grid to a source 
9j, j £ {1, • • • , fc}. It is easy to show that ||a;°||Q < fc, each 
column of A and B has unit norm, and f3° £ [— r, r]" with 

If-^. In addition, we let e = Rs with 



2n 



3m'J-10m^+7 



15 



such 



R = [Ri,--- ,Rk], and e = 
that ||e||2 < e (the information of fc and ||s||2 is used). The 
derivation for the setting of e is provided in Appendix F. Then 
the DOA estimation model can be written into the form of 
our studied model in (|9]l. The only differences are that A, 
B, x° and e are in the complex domain rather than the real 
domain and that e denotes a modeling error term rather than 
the measurement noise. It is noted that the robust stability 
results in SP-CS apply straightforward to such complex signal 
case with few modifications. The objective turns to recovering 
x° (its support actually) and /3°. According to Theorem [4] 
x" and (3° can be stably recovered if the D-RJP condition 
is satisfied. Denote x the recovered x", (3 the recovered /3°, 
and I the support of x. Then we obtain the recovered 9: 
9 — 9i + Kr^(3j^ where vi keeps only entries of a vector v 
on the index set I. The empirical results in Subsection |V-B| 
will illustrate the merits of applying the SP-CS framework to 
estimate DOAs. 
Remark 5: 
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(1) 



(2) 



Within the scope of DOA estimation, this work is related 
to spectral CS introduced in [^Fj. To obtain an accurate 
solution, the authors of |36| adopt a very dense sampling 
grid (that is necessary for any standard CS based methods 
according to the mentioned lower bound for the mean 
squared estimation error) and then prohibit a solution 
whose support contains near-located indices (that corre- 
spond to highly coherent columns in the overcomplete 
dictionary). In this paper we show that accurate DOA 
estimation is possible by using a coarse grid and jointly 
estimating the off-grid distance (the distance from a true 
DOA to its nearest grid point). 

The off-grid DOA estimation problem has been studied in 
fT2j , IpOl. The STLS solver in f20'l obtains a maximum a 
posteriori solution if (3° is Gaussian distributed. But such 
a condition does not hold in the off-grid DOA estimation 
problem. The SBI solver in fl^ proposed by the authors 
is based on the same model as in ([9]l and within the 



framework of Bayesian CS 1 37 1 



(3) 



The proposed P-BPDN can be extended to the multiple 
measurement vectors case like ^i-SVD to deal with DOA 
estimation with multiple snapshots. 

IV. Algorithms for P-BPDN 

A. Special Case: Positive Signals 

This subsection studies a special case where the original 
signal x° is positive-valued (except zero entries). Such a case 
has been studied in standard CS 138), |39|. By incorporating 



the positiveness of x°, P-BPDN is modified into the positive 
P-BPDN (PP-BPDN) problem 

||y-(A + BA)a;|| 

cc :>= 0, 

rl)^j3)p -rl, 



< e, 



min 1 X, subject to 



where is > with an elementwise operation and 0, 1 are 
column vectors composed of 0, 1 respectively with proper 
dimensions. It is noted that the robustly stable signal recovery 
results in the present paper apply directly to the solution to PP- 
BPDN in such case. This subsection shows that the nonconvex 
PP-BPDN problem can be transformed into a convex one and 
hence its optimal solution can be efficiently obtained. Denote 
p = (3 Q X. A new, convex problem (Pi) is introduced as 
follows. 



{Pi 



min l a;, subject to 



a; :>= 0, 

rx p —rx. 



Theorem 6: Problems PP-BPDN and (Pi ) are equivalent in 
the sense that, if (a;*, /3*) is an optimal solution to PP-BPDN, 
then there exists p* = 13* Qx* such that {x* ,p* ) is an optimal 
solution to (Pi), and that, if {x* ,p*) is an optimal solution to 

(Pi), then there exists /3* with /3* = | ^ ^^J*' ^ ^' 

' [ 0, otherwise 

such that {x* ,j3*) is an optimal solution to PP-BPDN. 

Proof: We only prove the first part of Theorem |6] using 

contradiction. The second part follows similarly. Suppose that 

[x* ,p*) with p* = /3* X* is not an optimal solution 



to (Pi). Then there exists [x' ,p') in the feasible domain 

X 



of (Pi) such that ||a;'||^ < 
p'j/x'j, if x'j > 0; 



0, 



otherwise 



^. Define f3' as /3j = 
It is easy to show that (a;',/3') is 



a feasible solution to PP-BPDN. By ||a;'||^ < \\x*\\-^ we 
conclude that {x* ,f3*) is not an optimal solution to PP-BPDN, 
which leads to contradiction. ■ 
Theorem |6] states that an optimal solution to PP-BPDN can 
be efficiently obtained by solving the convex problem (Pi). 



B. AA-P-BPDN: Alternating Algorithm for P-BPDN 



For general signals, P-BPDN in ( 1 1 1 is nonconvex. A simple 
method is to solve a series of BPDN problems with 



arg min 1 1 a; j 1 1 , subject to 



y-{A 



/3 



arg mm 

/3e[-r,r]'' 



BA^^^j X 

y-{A 



2 

BA) 



(20) 
(21) 



starting from /3'^*'^ — 0, where the superscript indicates the 
jth iteration and A^^^ = diag ( f3^^ A . Denote AA-P-BPDN 



the alternating algorithm defined by (20 1 and (21 1. To analyze 
AA-P-BPDN, we first present the following two lemmas. 
Lemma 1: For a matrix sequence composed of 



fat matrices, let T)^ = 
and V* ^{v: \\y-^*v 



0), 



< e 



2 < e} with e > 0. If ^ **, 
as i +00, then for any v E V* there exists a sequence 
{t;(j)}°^^ with e V^^\ j = 1, 2, • • • , such that v^^'' v, 
as j — > +00. 

Lemma [T] studies the variation of feasible domains , 
j = 1, 2, • • • , of a series of BPDN problems whose sensing 
matrices j = 1,2,- ••, converge to It states that 

the sequence of the feasible domains also converges to V* in 
the sense that for any point in 2?*, there exists a sequence 
of points, each of which belongs to one 2?-', that converges 
to the point. To prove Lemma [T| we first show that it holds 
for any interior point of V* by constructing such a sequence. 
Then we show that it also holds for a boundary point ofD* by 
that for any boundary point there exists a sequence of interior 
points of V* that converges to it. The detailed proof is given 
in Appendix C. 

Lemma 2: An optimal solution x* to the BPDN problem 
in (j4]) satisfies that x* — 0, if ||y||2 < e, or \\y — ^x*\\2 — e, 
otherwise. 

Proof: It is trivial for the case where ||y||2 < £■ Consider 
the other case where ||y||, > e. Note first that x* ^ 0. We 
use contradiction to show that the equality ||y — $a;*||2 = 
e holds. Suppose that \\y — ^x*]]^ < e. Introduce f {9) — 
\\y - e^x*\\^. Then /(O) > e, and /(I) < e. There exists Oq, 
< 6'o < 1, such that / [Oq) — e since / [0) is continuous on 
the interval [0, 1]. Hence, x' = 6*03;* is a feasible solution to 
BPDN in (jiji. We conclude that x* is not optimal by ||a;'||^ = 
6*0 II 33* 111 < II 3;* II 1^ which leads to contradiction. ■ 
Lemma [2] studies the location of an optimal solution to the 
BPDN problem. It states that the optimal solution locates at 
the origin if the origin is a feasible solution, or at the boundary 
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of the feasible domain otherwise. This can be easily observed 
from Fig.[T| Based on Lemmas [T] and [2] we have the following 
results for AA-P-BPDN. 

Theorem 7: Any accumulation point {x*,(3*) of the se- 
quence is a stationary point of AA-P- 
BPDN in the sense that"* 



13' 



= arg min 1 1 a? I L , subject to 

as 

\\y-{A + BA*)x\\^<e, 
= arg min \\y — {A + BA.) x*\\2 



(22) 
(23) 



with A* = diag (/3*). 

Theorem 8: An optimal solution {x*,(3*) to P-BPDN in 
( [TT] i is a stationary point of AA-P-BPDN. 

Theorem I?! studies the property of the solution 

produced by AA-P-BPDN. It shows that [x^i\l3^^"> 
bitrarily close to a stationary point of AA-P-BPDN as the 
iteration index j is large enough]^ Hence, the output of AA- 
P-BPDN can be considered as a stationary point provided that 
an appropriate termination criterion is set. Theorem |8] tells that 
an optimal solution to P-BPDN is a stationary point of AA- 
P-BPDN. So, it is possible for AA-P-BPDN to produce an 
optimal solution to P-BPDN. The proofs of Theorems [7] and 
|8] are provided in Appendix D and Appendix E respectively. 

Remark 6: During the revision of this paper, we have 
noted the following formulation of P-BPDN that possibly can 
provide an efficient approach to an optimal solution to P- 
BPDN. Let where Xi >p 0, a;_ >p and 



a; I x_ 



0. Then we have I a; I 



X- where 



applies elementwise. Denote p 
can be cast as follows: 



(3 Q X. A convex problem 



min 1^ {xj, 

X-|_,CC_ ,p 



subject to < 



x+ ^ 0, 
a; :>= 0, 

r {x^ + X- 





x+ 






A B] 


a;_ 




< e 




. P . 




2 




r{x+ 





(24) 



The above convex problem can be considered as a convex 
relaxation of P-BPDN since it can be shown (like that in 
Theorem |6]l that an optimal solution to P-BPDN can be 
obtained based on an optimal solution to the problem in 



(24 1 incorporated with an additional nonconvex constraint 
x+ a;_ = 0. An interesting phenomenon has been observed 
through numerical simulations that an optimal solution to the 
problem in (24i still satisfies the constraint a;_|_ a;_ = 0. 



Based on such an observation, an efficient approach to P- 
BPDN is to firstly solve ( [24) i, and then check whether its 
solution, denoted by (^x*^_,x*_,p*), satisfies x+ a;_ = 0. 

'^It is shown in the proof of Theorem ^ in Appendix D that the sequence 
is bounded. And it can be shown, for example, using 
contradiction, that for a bounded sequence {aj}'^-^, there exists an accumu- 
lation point of {a,j}°^i such that aj is arbitrarily close to it as j is large 
enough. 



If it does, then {x*,f3*) is an optimal solution to P-BPDN 

where x* = xl - x*_ and /3* = ( ^^^^^ ^ 

^ ■' [0, otherwise. 

Otherwise, we may turn to AA-P-BPDN again. But we note 

that it is still an open problem whether an optimal solution to 

(24 1 always satisfies the constraint a;-|_0a;_ = 0. In addition. 



the convex relaxation in (24i does not apply to the complex 



signal case as in DOA estimation studied in Subsection III-F 



C. Effectiveness of AA-P-BPDN 

As reported in the last subsection, it is possible for AA- 
P-BPDN to produce an optimal solution to P-BPDN. But it 
is not easy to check the optimality of the output of AA-P- 
BPDN because of the nonconvexity of P-BPDN. Instead, we 
study the effectiveness of AA-P-BPDN in solving P-BPDN in 
this subsection with the concept of effectiveness as defined in 
Subsection III-E By Corollary [T] a good signal recovery x of 
x° is not necessarily an optimal solution. It requires only that 
(^x,f3^, where f3 denotes the recovery of /3°, be a feasible 



solution to P-BPDN and that ||a;|L < \\x° 



holds. 



Hi — ,, Mi shown 
in the proof of Theorem jvj in Appendix D, that ^a;*^^), /S'-'-'^ 
for any j > 1 is a feasible solution to P-BPDN and that the se- 
quence { II 3'^"'' 11]^} is monotone decreasing and converges. 
So, the effectiveness of AA-P-BPDN in solving P-BPDN can 
be assessed via numerical simulations by checking whether 
1 1 3'"^'^ 11]^ < II 33° 111 holds with x"^"^ denoting the output of 
AA-P-BPDN. The effectiveness of AA-P-BPDN is verified in 
Subsection |V-A| via numerical simulations, where we observe 



that the inequality 
(over 3700 trials). 



\x 



AA\ 



< ||a;°||j^ holds in all experiments 



V. Numerical Simulations 

A. Verification of the Robust Stability 

This subsection demonstrates the robustly stable signal 
recovery results of SP-CS in the present paper, as well as the 
effectiveness of AA-P-BPDN in solving P-BPDN in ([TT]), via 



numerical simulations. AA-P-BPDN is implemented in Matlab 



with problems in ( 20 1 and ( 2 1 1 being solved using C VX |40) . 

^lllx'^'ll -llx'^^i'll I 

AA-P-BPDN is terminated as ^-^ , n ^ < 1 x IQ-^ 

or the maximum number of iterations, set to 200, is reached. 
PP-BPDN is also implemented in Matlab and solved by CVX. 

We first consider general signals. The sparse signal case is 
mainly studied. The variation of the signal recovery error is 
studied with respect to the noise level, perturbation level and 
number of measurements respectively. Besides AA-P-BPDN 
for P-BPDN in SP-CS, performances of three other approaches 
are also studied. The first one assumes that the perturbation is 
known a priori and recovers the original signal x° by solving, 
namely, the oracle (0-) BPDN problem 



a;||^, subject to \\y - {A + BA° 



< e. 



The O-BPDN approach produces the best recovery result of 
SP-CS within the scope of £1 minimization of CS since it ex- 
ploits the exact perturbation (oracle information). The second 
one corresponds to the robust signal recovery of perturbed CS 
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Fig. 2. Signal and perturbation recovery eiTors witli respect to the noise level 
e with parameter settings (n, m, k, r) = (200, 80, 10, 0.1). Both signal and 
l3° recovery errors of AA-P-BPDN for P-BPDN in SP-CS are proportional 
to e. 



Fig. 3. Signal and perturbation recovery errors with respect to the 
perturbation level in terms of r with parameter settings {n,m,k,e) = 
(200, 80, 10, 0.5). The error of AA-P-BPDN for P-BPDN in SP-CS slowly 
increases with the perturbation level and is quite close to that of the ideal 
case of O-BPDN for a moderate perturbation. 



as described in Subsection |II-C| and solves N-BPDN in (|7]) 
where eE.x" ~ ||-BA°a;°||2 is used though it is not available 
in practice. The last one refers to the other approach to SP- 
CS that seeks for the signal recovery by solving TPS-BPDN 
in ( fTSI as discussed in Subsection |III-E| 

The first experiment studies the signal recovery error with 
respect to the noise level. We set the signal length n = 200, 
sample size m — 80, sparsity level fc = 10 and perturbation 
parameter r — 0.1. The noise level e varies from 0.05 to 2 
with interval 0.05. For each combination of (n, m, k, r, e), the 
signal recovery error, as well as (3° recovery error (on the 
support of x°), is averaged over i? = 50 trials. In each trial, 
matrices A and B are generated from Gaussian distribution 
and each column of them has zero mean and unit norm after 
proper scaling. The sparse signal x° is composed of unit spikes 
with random signs and locations. Entries of /3° are uniformly 
distributed in [— r, r]. The noise e is zero mean Gaussian 
distributed and then scaled such that ||e||, = e. Using the same 
data, the four approaches, including O-BPDN, N-BPDN, TPS- 
BPDN and AA-P-BPDN for P-BPDN, are used to recover x° 
respectively in each trial. The simulation results are shown in 
Fig. [2] It can be seen that both signal and /3° recovery errors 
of AA-P-BPDN for P-BPDN in SP-CS are proportional to 
the noise, which is consistent with our robustly stable signal 
recovery result in the present paper The error of N-BPDN 
grows linearly with the noise but a large error still exhibits 
in the noise free case. Except the ideal case of O-BPDN, our 
proposed P-BPDN has the smallest error 

The second experiment studies the effect of the struc- 
tured perturbation. Experiment settings are the same as those 
in the first experiment except that we set (n, m, fc, e) = 
(200,80,10,0.5) and vary r e {0.05, 0.1, 1}. Fig. [3] 
presents our simulation results. A nearly constant error is 
obtained using O-BPDN in standard CS since the perturbation 
is assumed to be known in O-BPDN. The error of AA- 
P-BPDN for P-BPDN in SP-CS slowly increases with the 
perturbation level and is quite close to that of O-BPDN for 
a moderate perturbation. Such a behavior is consistent with 
our analysis. Besides, it can be observed that the error of N- 
BPDN grows linearly with the perturbation level. Again, our 
proposed P-BPDN has the smallest error except O-BPDN. 

The third experiment studies the variation of the recovery 
error with the number of measurements. We set (n, fc, r, e) — 




20 40 60 80 100 20 40 60 80 100 
Number of measurements Number of measurements 



Fig. 4. Signal and perturbation recovery eiTors with respect to the number of 
measurements with parameter settings (ra, fc, r, e) = (200, 10, 0.1, 0.2). AA- 
P-BPDN for P-BPDN in SP-CS has the best performance except the ideal 
case of O-BPDN. 



(200,10,0.1,0.2) and vai-y m e {30,35,- •• ,100}. Simula- 
tion results are presented in Fig. |4] Signal recovery errors of 
all four approaches decrease as the number of measurements 
increases. Again, it is observed that O-BPDN of the ideal case 
achieves the best result followed by our proposed P-BPDN. 
For example, to obtain the signal recovery error of 0.05, about 
55 measurements are needed for O-BPDN while the numbers 
are, respectively, 65 for AA-P-BPDN and 95 for TPS-BPDN. 
It is impossible for N-BPDN to achieve such a small error in 
our observation because of the existence of the perturbation. 

We next consider a compressible signal that is generated by 
taking a fixed sequence {2.8843 • with n = 200, 

randomly permuting it, and multiplying by a random sign 
sequence (the coefficient 2.8843 is chosen such that the 
compressible signal has the same £2 norm as the sparse signals 
in the previous experiments). It is sought to be recovered from 
m — 70 noisy measurements with e = 0.2 and r = 0.1. Give 
experiment results in one instance as an example. The signal 
recovery error of AA-P-BPDN for P-BPDN in SP-CS is about 
0.239, while errors of O-BPDN, N-BPDN and TPS-BPDN are 
about 0.234, 0.361 and 0.314 respectively. 

For the special positive signal case, an optimal solution to 
PP-BPDN can be efficiently obtained. An experiment result is 
shown in Fig. |5] where a sparse signal of length n — 200, 
composed of fc = 10 positive unit spikes, is exactly recovered 
from TO = 50 noise free measurements with r = 0.1 by solving 
(Pi)- 
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PP-BPDN, signal recovery error =1 .4065e-008 
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PP-BPDN, p° recovery error =2.8083e-009 
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Fig. 5. Exact recoveiy of a positive sparse signal from noise-free mea- 
surements with {m,n,k,r,e) = (200,50,10,0.1,0). PP-BPDN is solved 
by solving (Pi). (3° and its recovery are shown only on the support of x°. 
Black circles: original signal and (3°; red stars: recoveries. 



B. Empirical Results of DOA Estimation 

This subsection studies the empirical performance of the 
application of the studied SP-CS framework in DOA estima- 
tion. We consider the case of n = 90 and k — 2. Numer- 

< 



ical calculations show that the D-RIP condition S^k i"^) 

+ r^) + 1^ in Theorem is satisfied if m > 145. 
hough it ceases to be a "compressed" sensing problem in 
the case m > n, it still makes sense in SP-CS since there are 
2n variables to be estimated and hence the P-BPDN problem 
is still underdetermined as m < 2n. As noted in Subsection 
the D-RIP condition can be possibly relaxed using recent 
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techniques in standard CS, which may reduce the required m 
value. In addition, a RIP condition is a sufficient condition for 
guaranteed signal recovery accuracy while its conservativeness 
in standard CS has been studied in [Alj . We next choose a 
much smaller m = 30 (r « 0.302 in such a case) and show 
the empirical performance of the proposed SP-CS framework 
on such off-grid DOA estimation. 

The experimental setup is as follows. In each trial, the 
complex source signal s is generated with both entries having 
unit amplitude and random phases. 6i and 02 are generated 
uniformly from intervals ^] and i^] respectively 
(5.1° 7.7° apart in the DOA domain). "p-BPDN is solved 
using AA-P-BPDN whose settings are the same as those in 
Subsection V-A Our experimental results of the estimation 
error 6 for both sources are presented in Fig. |6] where 
1000 trials are used. It can be seen that P-BPDN performs 
well on the off-grid DOA estimation. All estimation errors 
he in the interval [— ^, ^] with most very close to zero. To 
achieve a possibly comparable mean squared estimation error, 
a grid of length at least n — 360 has to be used in standard 
CS based methods according to the lower bound mentioned in 
Subsection |III-F| An example of performance of SP-CS and 
standard CS on DOA estimation is shown in Fig. [7] where 
the two approaches share the same data set and n — 360 is 
set in standard CS. From the upper two sub-figures, it can be 
seen that SP-CS performs well on both source signal and (3° 
recoveries. From the lower left one, however, it can be seen 
that two nonzero entries are presented in the recovered signal 
around the location of each source when using standard CS. 
Such a phenomenon is much clearer in the last sub-figure. 



Fig. 6. Histogram of 6 estimation error for both sources using P-BPDN for 
SP-CS. Statistics including mean, variance and mean squared error (MSE) are 
shown. 
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Fig. 7. Performance comparison of SP-CS and standard CS (SCS) on DOA 
estimation. Upper left: signal recovery in SP-CS; upper right: /3° recovery 
in SP-CS (shown only on the signal support); lower left: signal recovery in 
standard CS; lower right: signal amplitude versus d (near the location of 
source 1) in SP-CS and standard CS. 



where it can be observed that a single peak exhibits at a place 
very close to the true location of source 1 using the proposed 
SP-CS framework while two peaks occurs at places further 
away from the true source in standard CS. 

VI. Conclusion 

This paper studied the CS problem in the presence of mea- 
surement noise and a structured matrix perturbation. A concept 
named as robust stability for signal recovery was introduced. It 
was shown that the robust stability can be achieved for a sparse 
signal by solving an li minimization problem P-BPDN under 
mild conditions. In the presence of measurement noise, the 
recovery error is at most proportional to the noise level and the 
recovery is exact in the special noise free case. A general result 
for compressible signals was also reported. An alternating 
algorithm named as AA-P-BPDN was proposed to solve the 
nonconvex P-BPDN problem, and numerical simulations were 
carried out, verifying our theoretical analysis. A practical 
application in DOA estimation was studied and satisfactory 
estimation results were obtained. 

The simulation results of DOA estimation suggest that the 
RIP condition for the robust stability is quite conservative in 
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practice. One future work is to relax such a condition. In our 
problem formulation, the signal x° and (3° that determines 
the matrix perturbation are jointly sparse. While this paper 
focuses on extracting the information that x° is sparse and 
that each entry of /3° lies in a bounded interval, such joint 
sparsity is not exploited. Inspired by the recent works on block 
and structured sparsity, e.g., | |42[ , | [43| , one future direction 
is to take into account the joint sparsity information in the 
signal recovery process to obtain possibly improved recovery 
performance. Our studied perturbed CS problem is related to 
the area of dictionary learning for sparse representation ]44| , 
where there is typically no a priori known structure in the 
overcomplete dictionary and a large number of observation 
vectors are important to make the learning process succeed. 
The studied problem in this paper can be considered as a 
dictionary learning problem but with a known structure in 
the dictionary, which leads to some similarity between our 
optimization approach and algorithms for dictionary learning, 
e.g., if-SVD 1^ and MOD ||45). Due to the known structure, 
it has been shown in this paper that a single observation vector 
is enough to learn the dictionary with guaranteed performance. 
Further relations deserve future studies. 



Appendix A 
Proof of Theorem[3] 

X 



Denote z — 
the problem in 



(BQx 



and similarly define z° and z*. Then 
|TO|l can be rewritten into 

a;||g , subject to y = i&z. (25) 



mm 

V/3e[-r,r]" 



Let 6k = 5k {^) hereafter for brevity. 

First note that x° is /c-sparse and z° is 2fc-D-sparse. 
Since {x*,f3*) is a solution to the problem in (25i, we 



have ||a;*||Q < l|3;°||o < k and, hence, z* is 2fc-D-sparse. 
By y = *2;° = ^z* we obtain *(z°-2;*) = and 
thus z° — z* = by S^k < 1 and the fact that z° — z* 
is 4fc-D-sparse. We complete the proof by observing that 

^ ^ [3° Q x" - (3* Q X* 



0. 



Appendix B 
Proofs of Theorems|4]and[5] 

We only present the proof of Theorem |5] since Theorem |4] 
is a special case of Theorem [5] We first show the following 
lemma. 

Lemma 3: We have 



And thus 



|(*i;,*t;')| < 



< 6. 



2{k+k'), 



which completes the proof. ■ 
Using the notations z, z°, z* and 5k in Appendix A, P- 
BPDN in ( [TT| i can be rewritten into 



mm 



\x\\^ , subject to \\y ~ ^z\\^ < e. (26) 



Let h ~ X* — x° and decompose h into a sum of fc-sparse 
vectors /ijo , , , • • s where Tq denotes the set of indices 
of the k largest entries (in absolute value) of x°, Ti the set of 
the k largest entries of hr^ with Tq being the complementary 
set of Tq, T2 the set of the next k largest entries of hxg and so 



on. We abuse notations z^, = 



(3* a;* 



J =0,1,2,. 



and similarly define z^^ . Let f — z* —z° and jj,, — zip, —z"^. 



for j = 0,1,2, 



For brevity we write Tq 



To U Ti. 



To bound \\h\\^, in the first step we show that ||/i'T,5:J|2 is 
essentially bounded by ||/itoiII2' ^'^'^ '■'^^'^ ™ second step 
we show that ||ft.Toill2 sufficiently small. 

The first step follows from the proof of Theorem 1.3 in Q. 
Note that 

\\hT,\\^<k^'^\\hT^\\^<k-^'^\\hT,_,\\^, J>2, (27) 
and thus 



I -"oi 112 



J>2 



< 



E 

J>2 



3 II2 



(28) 



Since x* = x° + h is an optimal solution, we have 



c°||i >\\x° + h\\,= J2 kj+'^jl 



> X°n 



E 



and thus 



To 111 



(29) 



(30) 



\{^v,^v')\<52ik+k')\\v\\Jv% 

for all 2fc-D-sparse v and 2A:'-D-sparse v' supported on 
disjoint subsets. 

Proof: Without loss of generality, assume that v and v' 
are unit vectors with disjoint supports as above. Then by the 
definition of D-RIP and \\v ± v'\\l = \\v\\l + \\v'\\l = 2 we 
have 

2 (1 - S^ik+k')) <\\^v± ^v'Wl < 2 (1 + 52(fe+fc')) • 



By (|28|, (|30| and the inequality ||/itoIIi ^ ^^^^ W^Toh we 
have 



I -'oi II2 — / ' 



\hT,\\^ < \\hT„\\2 + ' 



2k~^'^eo (31) 

with eo = ||a;° — x'^^^. 

In the second step, we bound ||/itoiII2 utilizing its 
relationship with \\fToi \\2' Note that /j.^ for each j = 0, 1, • • • 
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is 2fc-D-spai-se. By */toi = */ - J2j>2 ^fr, we have 

< K*/to„*/)|+E|(*/to,*/t, 

< ll*/ToJl2-|l*/ll2 + ^"4/.||/ToLE|kT, 



By ([38]l and the RIP we have 

<I|S[(/3to-/3to)0^To]L 
r n 1 



j>2 



+^4k ii/tiLE k^. 



J>2 



< ||/3.„J|2 2eVl + ^4fc + \/2^4fcE|kT, 



J>2 



(32) 



(33) 



(/3to 

<II*/Il2 + 

<2e+ VT + 



'To 



a;° 



h 

2||*IUII^ 



B 



13; 



/3yc ^ X^a 



2r||B|Leo 



<C4e + (^c^k + cg^ eo 



We used Lemma 3] in (f32|. In ( |33| ), we used the D-RIP, and 
inequalities H/toIIz + II/tiL ^ ^ll-^ToiL 

- II* {z* - z")\\^ <\\y- ^z*\\^ + \\y - ^z% < 2e. 



with Sk (B) 

2Vl+'-2||*|l2Ci 
1-cn 



< 

2C2 



1-Co 



2r ) 11*11^, and thus 



^4fe 

+ 2 



C4 
1*11 



and cg = 



By noting that /3°,/3* e [-r, r]" and 



/3* + (/3* - 13°) a;' 



(34) 
(35) 



<- 



II(/3to 
1 



bll2 



34fe 



C4e + ^csfc 



-1/2 



+ Cg eo 



we have 



which concludes ( [T5| ). We complete the proof by noting that 
the above results make sense if cq < 1, i.e., 

1 



< y^l + r^\\hT^\\^ + 2r 



Meanwhile, 



E \rh , - E \\^% 



i>2 



i>2 



Xn 



(36) 

< eo. (37) 



1 



Applying the D-RIP, (|36jl and then (|3T]l and (|37ji it gives 

|2 



(l-<54fe)||/Toj2<r/T„,||2 



<l|/ 



Tai II 2 



{2\/l + 54fce+ V2(l + r2)^4fc ||/itoII 



Appendix C 
Proof of Lemma[T] 

We first consider the case where v is an interior point of 
V*, i.e., it holds that ||y — **i'||2 = eo < £■ Let r/ = e — eq. 
Construct a sequence {v'^-'^} such that \\v^^'> — v\\^ < 
It is obvious that v^^^^ — !• v. We next show that t;'^^ e as 
j is large enough. By t;'^^ ^ v and that the 

sequence {v'-'^} is bounded, there exists a positive integer 
jo such that, as j > jo. 



and thus 



||/lToJl2 < II/T0JI2 
<Cie + Co||/lToJl2+ fc2fc"^/^ 



II** 



,(i) 



t> — V 



< W2 



with Co 



^2(l+r^)g4fc 



1-54 



, Ci 



C3 eo 



2^1 + ^4 

1— 541, 



-, C2 = 2co and C3 



272^4)='- ^ Hence, we get a bound 

\\hT0A\2 < (l-co)^^ cie+ (c2fc-i/2_^C3 
which together with pT[ ) gives 

||/l|l2< l|/lToi 112+11^1 



$0') I I, 



(i) 



* 



eo 



** 



** 
* 



01 II2 



<2||/lToJl2 

2ci 

1 - Co 



2/fc'i/2 



2C2 



eo 



1 - Co 



2U-1/2 



2C3 



1-co 



(38) 



eo, 



which concludes ( [T4] i. 



Hence, as j > jo^ 

(y - **t)) -+ 

< \\y - ^*v\\^ 4 

< ||y - **'i;||2 4 
<eo + 77/2 + 77/2 = e, 

from which we have v*^^' G I?-' for j > jo- By re-selecting 
arbitrary v^^^ e I?^ for j < jo we obtain the conclusion. 

For the other case where 1; is a boundary point of 2?*, there 
exists a sequence C V* with all V(i-j being interior 

points of V* such that v, as I ^ +00. According to 



V — v^-'^ 
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sequence 



the first part of the proof, for each I = 1, 2, • • • , there exists a 
^v'-fj^ I with v'-fj^^ e I?-', j = 1, 2, • • • , such that 

<^j,jj —r as j +00. The sequence is what 

we expected since 



Appendix E 
Proof of Theorem[8] 

We need to show that an optimal solution {x*,f3*) satisfies 
221) and (|23ll. It is obvious for (|22li. For (|23ll, we discuss 



as J 



U) 

+00. 



< 



'U) 



'U) 



0, 



Appendix D 
Proof of Theorem|7] 

We first show the existence of an accumulation point. It 
follows from the inequality 



two cases based on Lemma |2] If |jy||2 < e, then x* = 
and, hence, (23 1 holds for any (3* S [— r, r]". If ||y|L > e, 
\\y - {A + BA*) x*\\2 = e holds by ([22] i and Lemma|2| Next 
we use contradiction to show that (|23 i holds in such case. 

Suppose that (23 i does not hold as jyHj > e. That is, there 
exists (3' E [— r, such that 



\y-(A + BA') 



X 



< \\y-{A + BA*)x*\\^ = e 



holds with A' = diag (/3'). Then by Lemma|2]we see that x* 
is a feasible but not optimal solution to the problem 



y - (a + BA'-J^ 



.0) 



< 



y (a + BA^^-^A x^^^ 



< e 



min||a;||-^, subject to (A 



BA') x\\.^ < e. 



(40) 



that a;^^^ is a feasible solution to the problem in (20i, and 
thus ||a;'^-'+-'-^||j^ < ||a;*^-'^||^ for j = 1,2.---. Then we have 



.(i)| 



A^y 



for j = 1, 2, • • • , since A^y is 



a feasible solution to the problem in (|20]) at the first iteration 
with the superscript ^ denoting the pseudo-inverse operator 
This together with /S*-^-' G [— r, r]", j — 1,2,---, leads to 

that the sequence | ^a;^^\ /3'-'^^ | is bounded. Thus, there 

exists an accumulation point (a;*,/3*) of | |^a;^^\ /j'-'-'^ | 

For the accumulation point (a;*,/3*) there exists a subse- 
quence |(a;(J'),^(^'))|^^ of \^[x''i\l3^^^^Y° ^ such that 

all /3 e [-r, r]", 

< y - [A + BA)x^^''i 
at both sides of which by taking / — > +oo, we have, for all 

\\y-{A + BA*) x*\\^ < \\y -{A + BA) x*\\, , 
which concludes (|2T 



as I 



-oo. By (21 



we have, for 



y (a + BA'^^'A x^^'^ 



For (22 1, we first point out that 



as J 



-oo, since |||a;(-''|L| is decreasing and 

L II 111 J j = l ° 

X* is one of its accumulation points. As in Lemma u\ 



let 



|a; : y - (^A + B A'-^^^ x < e| 



and V* 



y - \^A + BA^^' J X 

{x : \\y - (A + BA*)x\\^ < e}. By A + BA^^' 
BA*, as I +00, and Lemma[T] for any x £ V* there exists 
a sequence with a;(;') G P^', / = 1, 2, • • • , such that 

Xf^i-f X, as I ^ +00. By (j20|l, we have, for Z = 1, 2, • • • , 



< \\x 



at both sides of which by taking I 



(Olli ' 

oo, we have 



X 



li < lla^lli 



(39) 



since ||a;(^^||j^ \\x*\\j^, as j +oo, and a;^;) x, as 
I +00. Finally, (22i is concluded as (39i holds for arbitrary 

X e V*. 



Hence, ||a;'||j< ||a;*||]^ holds for an optimal solution x' to the 
problem in (40 1. Meanwhile, (a;',/3') is a feasible solution to 



the P-BPDN problem in ( 11 1. Thus (x* , f3*) is not an optimal 
solution to the P-BPDN problem in (flTi by ||a;'||^ < ||a;*||p 
which leads to contradiction. 

Appendix F 
Deriavation of e in Subsection IV-BI 

By ( 19 1, we have for ^ = 1, • • • ,m, j — 1, ■ ■ ■ , fc. 



Rij - 



1j - Oi, 



(41) 



where ^ is between 9j and 6i., A'{^ (^) 
Thus, we have for j = 1, • • • , fc. 



< 



15 



Finally, it gives the expression of e by observing that 

||e|U-||fis|U<ll«llFll^ll.<^ll^i|l2 11^112- 



(42) 



(43) 
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